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Abstract 

We present a new algorithm for Independent Component Analysis (ICA) which has provable per- 
formance guarantees. In particular, suppose we are given samples of the form y = Ax + r\ where A is 
an unknown n x n matrix and i is a random variable whose components are independent and have a 
fourth moment strictly less than that of a standard Gaussian random variable and n is an n-dimensional 
Gaussian random variable with unknown covariance E: We give an algorithm that provable recovers A 
and E up to an additive e and whose running time and sample complexity are polynomial in n and 1/e. 
To accomplish this, we introduce a novel "quasi-whitening" step that may be useful in other contexts 
in which the covariance of Gaussian noise is not known in advance. We also give a general framework 
for finding all local optima of a function (given an oracle for approximately finding just one) and this 
is a crucial step in our algorithm, one that has been overlooked in previous attempts, and allows us to 
control the accumulation of error when we find the columns of A one by one via local search. 

1 Introduction 

We present an algorithm (with rigorous performance guarantees) for a basic statistical problem. Suppose 
X] is an independent n-dimensional Gaussian random variable with an unknown covariance matrix S and A 
is an unknown n x n matrix. We are given samples of the form y = Ax + r\ where a; is a random variable 
whose components are independent and have a fourth moment strictly less than that of a standard Gaussian 
random variable. The most natural case is when x is chosen uniformly at random from {+1, —1}™, although 
our algorithms in even the more general case above. Our goal is to reconstruct an additive approximation to 
the matrix A and the covariance matrix E running in time and using a number of samples that is polynomial 
in n and -, where e is the target precision (see Theorem 1.1) This problem arises in several research directions 
within machine learning: Independent Component Analysis (ICA), Deep Learning, Gaussian Mixture Models 
(GMM), etc. We describe these connections next, and known results (focusing on algorithms with provable 
performance guarantees, since that is our goal). 

Most obviously, the above problem can be seen as an instance of Independent Component Analysis (ICA) 
with unknown Gaussian noise. ICA has an illustrious history with applications ranging from econometrics, 
to signal processing, to image segmentation. The goal generally involves finding a linear transformation of 
the data so that the coordinates are as independent as possible [1, 2, 3]. This is often accomplished by finding 
directions in which the projection is "non-Gaussian" [4]. Clearly, if the datapoint y is generated as Ax (i.e., 
with no noise r\ added) then applying linear transformation A^ 1 to the data results in samples A~ 1 y whose 
coordinates are independent. This restricted case was considered by Comon [1] and Frieze, Jerrum and 
Kannan [5] , and their goal was to recover an additive approximation to A efficiently and using a polynomial 
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number of samples. (We will later note a gap in their reasoning, albeit fixable by our methods. See also 
recent papers by Anandkumar et al., Hsu and Kakade[G, 7], that do not use local search and avoids this 
issue.) To the best of our knowledge, there arc currently no known algorithms with provable guarantees for 
the more general case of ICA with Gaussian noise (this is especially true if the covariance matrix is unknown, 
as in our problem), although many empirical approaches are known, (eg. [8], the issue of "empirical" vs 
"rigorous" is elaborated upon after Theorem 1.1.) 

The second view of our problem is as a concisely described Gaussian Mixture Model. Our data is generated 
as a mixture of 2™ identical Gaussian components (with an unknown covariance matrix) whose centers are 
the points {Ax : x G { — 1,1}™}, and all mixing weights are equal. Notice, this mixture of 2™ Gaussians 
admits a concise description using 0(n 2 ) parameters. The problem of learning Gaussian mixtures has a long 
history, and the popular approach in practice is to use the EM algorithm [9], though it has no worst-case 
guarantees (the method may take a very long time to converge, and worse, may not always converge to the 
correct solution). An influential paper of Dasgupta [10] initiated the program of designing algorithms with 
provable guarantees, which was improved in a sequence of papers [11, 12, 13, 14]. But in the current setting, 
it is unclear how to apply any of the above algorithms (including EM) since the trivial application would 
keep track of exponentially many parameters - one for each component. Thus, new ideas seem necessary to 
achieve polynomial running time. 

The third view of our problem is as a simple form of autoencoding [15]. This is a central notion in Deep 
Learning, where the goal is to obtain a compact representation of a target distribution using a multilayered 
architecture, where a complicated function (the target) can be built up by composing layers of a simple 
function (called the autoencodcr [16]). The main tenet is that there are interesting functions which can be 
represented concisely using many layers, but would need a very large representation if a "shallow" architecture 
is used instead). This is most useful for functions that are "highly varying" (i.e. cannot be compactly 
described by piecewise linear functions or other "simple" local representations). Formally, it is possible to 
represent using just (say) n 2 parameters, some distributions with 2" "varying parts" or "interesting regions." 
The Restricted Boltzmann Machine (RBM) is an especially popular autoencoder in Deep Learning, though 
many others have been proposed. However, to the best of our knowledge, there has been no successful 
attempt to give a rigorous analysis of Deep Learning. Concretely, if the data is indeed generated using 
the distribution represented by an RBM, then do the popular algorithms for Deep Learning [17] learn the 
model parameters correctly and in polynomial time? Clearly, if the running time were actually found to 
be exponential in the number of parameters, then this would erode some of the advantages of the compact 
representation. 

How is Deep Learning related to our problem? As noted by Freund and Haussler [18] many years ago, 
an RBM with real-valued visible units (the version that seems more amenable to theoretical analysis) is 
precisely a mixture of exponentially many standard Gaussians. It is parametrized by annxm matrix A 
and a vector 9 £ R™. It encodes a mixture of n-dimensional standard Gaussians centered at the points 
{Ax : x € { — 1, 1} }j where the mixing weight of the Gaussian centered at Ax is exp(||Ax||| + 6 ■ x). This 
is of course reminiscent of our problem. Formally, our algorithm can be seen as a nonlinear autoencoding 
scheme analogous to an RBM but with uniform mixing weights. Interestingly, the algorithm that we present 
here looks nothing like the approaches favored traditionally in Deep Learning, and may provide an interesting 
new perspective. 

1.1 Our results and techniques 

We give a provable algorithm for ICA with unknown Gaussian noise. We have not made an attempt to 
optimize the quoted running time of this model, but we emphasize that this is in fact the first algorithm 
with provable guarantees for this problem and moreover we believe that in practice our algorithm will run 
almost as fast as the usual ICA algorithms, which arc its close relatives. 

Theorem 1.1 (Main, Informally). There is an algorithm that recovers the unknown A and E up to additive 
error e in each entry in time that is polynomial in n, || A||2, ||£||2; V £ > V^mmC-A) where || • H2 denotes the operator 
norm and A m j n (-) denotes the smallest eigenvalue. 
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The classical approach for ICA initiated in Comon [1] and Frieze, Jerrum and Kannan [5]) is for the 
noiseless case in which y = Ax. The first step is whitening, which applies a suitable linear transformation 
that makes the variance the same in all directions, thus reducing to the case where A is a rotation matrix. 
Given samples y = Rx where R is a rotation matrix, the rows of R can be found in principle by computing 
the vectors u that are local minima of E\(u ■ y) 4 ]. Subsequently, a number of works (see e.g. [19, 20]) 
have focused on giving algorithms that arc robust to noise. A popular approach is to use the fourth order 
cumulant (as an alternative to the fourth order moment) as a method for "denoising," or any one of a number 
of other functionals whose local optima reveal interesting directions. However, theoretical guarantees of these 
algorithms are not well understood. 

The above procedures in the noise- free model can almost be made rigorous (i.e., provably polynomial 
running time and number of samples), except for one subtlety: it is unclear how to use local search to find all 
optima in polynomial time. In practice, one finds a single local optimum, projects to the subspace orthogonal 
to it and continues recursively on a lower-dimensional problem. However, a naive implementation of this idea 
is unstable since approximation errors can accumulate badly, and to the best of our knowledge no rigorous 
analysis has been given prior to our work. (This is not a technicality: in some similar settings the errors are 
known to blow up exponentially [21].) One of our contributions is a modified local search that avoids this 
potential instability and finds all local optima in this setting. (Section 4.2.) 

Our major new contribution however is dealing with noise that is an unknown Gaussian. This is an 
important generalization, since many methods used in ICA are quite unstable to noise (and a wrong estimate 
for the covariance could lead to bad results). Here, we no longer need to assume we know even rough estimates 
for the covariance. Moreover, in the context of Gaussian Mixture Models this generalization corresponds to 
learning a mixture of many Gaussians where the covariance of the components is not known in advance. 

We design new tools for denoising and especially whitening in this setting. Denoising uses the fourth 
order cumulant instead of the fourth moment used in [5] and whitening involves a novel use of the Hessian 
of the cumulant. Even then, we cannot reduce to the simple case y = Rx as above, and are left with a 
more complicated functional form (see "quasi-whitening" in Section 2.) Nevertheless, we can reduce to an 
optimization problem that can be solved via local search, and which remains amenable to a rigorous analysis. 
The results of the local optimization step can be then used to simplify the complicated functional form and 
recover A as well as the noise E. We defer many of our proofs to the supplementary material section, due 
to space constraints. 

In order to avoid cluttered notation, we have focused on the case in which x is chosen uniformly at 
random from { — 1, +1}™, although our algorithm and analysis work under the more general conditions that 
the coordinates of x are (i) independent and (ii) have a fourth moment that is less than three (the fourth 
moment of a Gaussian random variable). In this case, the functional P(u) (see Lemma 2.2) will take the 
same form but with weights depending on the exact value of the fourth moment for each coordinate. Since 
we already carry through an unknown diagonal matrix D throughout our analysis, this generalization only 
changes the entries on the diagonal and the same algorithm and proof apply. 

2 Denoising and quasi-whitening 

As mentioned, our approach is based on the fourth order cumulant. The cumulants of a random variable 
are the coefficients of the Taylor expansion of the logarithm of the characteristic function [22]. Let n r (X) 
be the r th cumulant of a random variable X. We make use of: 

Fact 2.1. (i) If X has mean zero, then K±{X) = E[A 4 ] — 3E[A 2 ] 2 . (ii) If X is Gaussian with mean \i and 
variance a 2 , then K\(X) = fx, k 2 (A) = a 2 and K r (X) = for all r > 2. (Hi) If X and Y are independent, 
then n r (X +Y) = K r (X) + n r (Y). 

The crux of our technique is to look at the following functional, where y is the random variable Ax + n 
whose samples are given to us. Let u € M. n be any vector. Then P(u) = —K±(u T y). Note that for any 
u we can compute P(u) reasonably accurately by drawing sufficient number of samples of y and taking 
an empirical average. Furthermore, since x and rj are independent, and r\ is Gaussian, the next lemma is 
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immediate. We call it "denoising" since it allows us empirical access to some information about A that is 
uncorrupted by the noise rj. 

Lemma 2.2 (Denoising Lemma). P(u) = 2 ^" =1 (u T A)^. 

Proof: The crucial observation is that u T y = u T Ax + u t t] is the sum of two independent random variables, 
Ax and rj and that P(u) = —kh(u t Ax + u T r}) = — k,^{u t Ax) — Ki(u T r}) = —k^(u t Ax). So in fact, the 
functional P{u) is invariant under additive Gaussian noise independent of the variance matrix E. This 
vastly simplifies our computation: 

n 

E[(u T Axf] = ^(u T A)i Efctf + 6 Y,(u T A)l(u T A)> E[xf] E^] 

i—1 i<j 
n n 

= J> T A)J + 6^^^^ = ~2]>> T ^ + Z{u T AA T uf 

i—1 i<Cj i—1 

Furthermore E[(w T ylx) 2 ] 2 = (u T AA T u) 2 and we conclude that 

n 

P(u) = -K 4 (u T y) = -E[(/ii) 4 ]+3E[(/Aa;) 2 ] 2 = 2^(u T A)f. 

i=l 



2.1 Quasi-whitening via the Hessian of P(u) 

In prior works on ICA, whitening refers to reducing to the case where y = Rx for some some rotation matrix 
R. Here we give a technique to reduce to the case where y = RDx + rf where rj' is some other Gaussian 
noise (still unknown), R is a rotation matrix and D is a diagonal matrix that depends upon A. We call this 
quasi-whitening. Quasi- whitening suffices for us since local search using the objective function K4(u T y) will 
give us (approximations to) the rows of RD, from which we will be able to recover A. 

Quasi- whitening involves computing the Hessian of P{u), which recall is the matrix of all 2nd order 
partial derivatives of P(u). Throughout this section, we will denote the Hessian operator by JC. In matrix 
form, the Hessian of P(u) is 

n2 n n 

P(u)^2A^2A :k A J , k (A k -uj 2 ; n(P{U))^2AY^{A k -u) 2 A k A T k ^AD A {u)A T 



1 ■> fe=i fe=i 

where A k is the fc-th column of the matrix A (we use subscripts to denote the columns of matrices throught 
the paper). Da(u) is the following diagonal matrix: 

Definition 2.3. Let Da(u) be a diagonal matrix in which the k th entry is 2A(A k ■ u) 2 . 

Of course, the exact Hessian of P{u) is unavailable and we will instead compute an empirical approxima- 
tion P(u) to P(u) (given many samples from the distribution), and we will show that the Hessian of P(u) 
is a good approximation to the Hessian of P{u). 

Definition 2.4. Given 2N samples y\, y' x , t/2, y'i— , VN, Vn 01 the random variable y, let 

p M - w I> T ^ 4 + 1 f> T y*) 2 (^) 2 - 

i=l i=l 

Our first step is to show that the expectation of the Hessian of P{u) is exactly the Hessian of P(u). In 
fact, since the expectation of P(u) is exactly P(u) (and since P(u) is an analytic function of the samples 
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and of the vector u), we can interchange the Hessian operator and the expectation operator. Roughly, one 
can imagine the expectation operator as an integral over the possible values of the random samples, and as 
is well-known in analysis, one can differentiate under the integral provided that all functions are suitably 
smooth over the domain of integration. 

Claim 2.5. E y y [-(u T y) 4 + 3(u T y) 2 (u T y') 2 ] = P{u) 

This claim follows immediately from the definition of P(u), and since y and y' are independent. 
Lemma 2.6. M(P(u)) = V y . y ,[K{-{u T yf + 'S(u T y) 2 (u T y') 2 )} 

Next, we compute the two terms inside the expectation: 
Claim 2.7. 3i((u T y) 4 ) = \2(u T y) 2 yy T 

Claim 2.8. % {{u T y) 2 {u T y') 2 ) = 2(u T y') 2 yy T + 2(u T y) 2 y> \y') T + A{u T y){u T y'){y{y') T + {y')y T ) 

Let X m i n {A) denote the smallest eigenvalue of A. Our analysis also requires bounds on the entries of 
D A (u ): 

Claim 2.9. If uq is chosen uniformly at random then with high probability for all i, 

min IIAJInTT. -4 < Da(uo)h < maxIMJIn — - — 
j=i ' i=i n 

Proof: We can bound max™ =1 \Aj ■ u\ by max™ =1 || A t || 2 thus the bound for m&x™ =1 (DA(u )) it i follows. 
Note that with high probability the minimum absolute value of n Gaussian random variables is at least 1/n 2 , 
hence min™ =1 (.D,i(uo))i,i > mm " =1 j| Aijl^n -4 . ' 

Lemma 2.10. IfuQ is chosen uniformly at random and furthermore we are given 2N = poly(rt, 1/e, l/X m in{A), 
|| A\\2, ||S||2) samples of y, then with high probability we will have that (1 — £)ADa(uq)A t ^ "K(P{uq)) < 
(l + e)AD A (u )A T . 

Proof: First we consider each entry of the matrix updates. For example, the variance of any entry in 
Jf((w T y) 4 ) = \2(u T y) 2 yy T can be bounded by , which we can bound by E[||y|||] < 0(E[|| Aa;||| + WvWl])- 
This can be bounded by 0(n 4 (|| ^4||| + ||Sj||)). This is also an upper bound for the variance (of any entry) 
of any of the other matrix updates when computing J£(P(uo)). 

Applying standard concentration bounds, poly(n, 1/e', ||^4||2, ll^lb) samples suffice to guarantee that all 
entries of JC(P(tio)) are e' close to JC(P(u)). The smallest eigenvalue of J£(P(u)) = ADa(uo)A t is at least 
Xmin(A) 2 min™ =1 j| Ai\\\n~ A where here we have used Claim 2.9. If we choose e' = poly(l/n, X m in{A), e), then 
we are also guaranteed (1 - e)ADA{u )A T ^ Jf(P(u )) ^ (1 + e)ADA{u )A T holds. ■ 

Lemma 2.11. Suppose that (1 - e)AD A (u )A T X M ^ (1 + e)AD A (u )A T , and let M = BB T . Then there 
is a rotation matrix R* such that \\B~ 1 ADa{uo) 1 ^ 2 — R*\\f < \fnt. 

The intuition is: if any of the singular values of B^ 1 ADa(uq) 1 ^ 2 are outside the range [1 — e, 1 + e], we 
can find a unit vector x where the quadratic forms x t ADa{uo)A t x and x T Mx are too far apart (which 
contradicts the condition of the lemma). Hence the singular values of B~ 1 ADa(uo) 1 ^ 2 can all be set to one 
without changing the Froebenius norm of B~ 1 ADa(uq) 1 ^ 2 too much, and this yields a rotation matrix. 
Proof: Let M = AD A {u )A T and let C = AD A {u ) 1/2 , and so M = CC T and M = BB T . The condition 
(1— e)M < M -< (l+e)M is well-known to be equivalent to the condition that for all vectors x, (l — e)x T Mx < 
x T Mx < (1 + e)x T Mx. 

Suppose for the sake of contradiction that S = B~ 1 C has a singular value outside the range [1 — e, 1 + 
e]. Assume (without loss of generality) that S has a singular value strictly larger than 1 + e (and the 
complementary case can be handled analogously). Hence there is a unit vector y such that y T SS T y > 1 + e. 
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But since BSS T B T = CC T , if we set x T = y T B^ 1 then we have x T Mx = x T BB T x = y T y = 1 but 
x T Mx = x T CC T x = x T BSS T B T x = y T SS T y > 1 + e. This is a contradiction and so we conclude that all 
of the singular values of B~ 1 C are in the range [1 — e, 1 + e]. 



Let WEV T be the singular value decomposition of B~ X C. If we set all of the diagonal entries in S to 1 we 
obtain a rotation matrix R* = UV T . And since the singular values of B~ 1 C arc all in the range [1 — e, 1 + e], 
we can bound the Froebenius norm of B~ X C — R*: \\B~ 1 C — R*\\f < \fne, as desired. ■ 

3 Our algorithm (and notation) 

In this section we describe our overall algorithm. It uses as a blackbox the denoising and quasi-whitening 
already described above, as well as a routine for computing all local maxima of some "well-behaved" functions 
which is described later in Section 4. 

Notation: Placing a hat over a function corresponds to an empirical approximation that we obtain from 
random samples. This approximation introduces error, which we will keep track of. 

Step 1: Pick a random uq € K" and estimate the Hessian ^K(P(uq)). Compute B such that J-C(P(uo)) = 
BB T . Let D = Da{uq) be the diagonal matrix defined in Definition 2.3. 
Step 2: Take 2N samples yi,y2, ■■■,yN,y'i,y2, ■■■,y / N> andlet 



which is an empirical estimation of P'(u). 

Step 3: Use the procedure AllOPT(P'(u), /3, 5', /?', 5') of Section 4 to compute all n local maxima of the 
function P'{u). 

Step 4: Let R be the matrix whose rows are the n local optima recovered in the previous step. Use procedure 
Recover of Section 5 to find A and S. 

Explanation: Step 1 uses the transformation B~ x computed in the previous Section to quasi- whiten the 
data. Namely, we consider the sequence of samples z — B~~ 1 y, which are therefore of the form R'Dx + rj' 
where r\ = B~ x r\, D = Da(uq) and R' is close to a rotation matrix R* (by Lemma 2.11). In Step 2 we 
look at K4((u T z)), which effectively denoises the new samples (see Lemma 2.2), and thus is the same as 
KiiR'D^^x). Let P'(u) = k 4 {u t z) = K 4 {u T B~ 1 y) which is easily seen to be E[{u T R' D~ 1 / 2 x) i }. Step 2 
estimates this function, obtaining P'{u). Then Step 3 tries to find local optima via local search. Ideally we 
would have liked access to the functional P*(u) = (u T R*x) 4 since the procedure for local optima works only 
for true rotations. But since R' and R* arc close we can make it work approximately with P'(u), and then 
in Step 4 use these local optima to finally recover A. 

Theorem 3.1. Suppose we are given samples of the form y = Ax + r] where x is uniform on {+1, —1}™, A is 
an nx n matrix, rj is an n-dimensional Gaussian random variable independent of x with unknown covariance 
matrix S. There is an algorithm that with high probability recovers \\A — AHdiag(ki)\\p < e where H is some 
permutation matrix and each ki G {+1,-1} and also recovers ||S — S||f < £■ Furthermore the running time 
and number of samples needed are poly(n, 1/e, ||A|| 2 , ||S|| 2 , l/A m ,; n (A)) 

Proof: In Step 1, by Lemma 2.11 we know once we use z = B _1 y, the whitened function P'{u) is inverse 
polynomially close to P*(u). Then by Lemma 5.3, the function P'(u) we get in Step 2 is inverse polynomially 
close to P'(u) and P*(u). Theorem 4.6 and Lemma 5.5 show that given P'(u) inverse polynomially close 
to P*(u), Algorithm 2: : AllOPT finds all local maxima with inverse polynomial precision. Finally by 
Theorem 5.6 we know A and W are recovered correctly up to additive e error in Frobcnius norm. The running 
time and sampling complexity of the algorithm is polynomial because all parameters in these Lemmas are 
polynomially related. H 

Note that here we recover A up to a permutation of the columns and sign-flips. In general, this is all we 
can hope for since the distribution of x is also invariant under these same operations. Also, the dependence 
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Algorithm 1. LocalOPT, Input:/(u), u s , (3, S Output: vector v 



1. Set u -<— u s - 

2. Maximize (via Lagrangian methods) Proj ±u (V/(«)) T £ + ^ T Proj Xu (3i(f(u)))Z - \ (£/(«)) ■ 
Subject to ||£|| 2 < /3' and u T £ = 

3. Let £ be the solution, u — ic^|ti 

4. If /(ft) > /(w) + 5/2, set it «- u and Repeat Step 2 

5. Else return it 



of our algorithm on the various norms (of A and E) seems inherent since our goal is to recover an additive 
approximation, and as we scale up A and/or E, this goal becomes a stronger relative guarantee on the error. 

4 Framework for iteratively finding all local maxima 

In this section, we first describe a fairly standard procedure (based upon Newton's method) for finding 
a single local maximum of a function /* : W — > M among all unit vectors and an analysis of its rate of 
convergence. Such a procedure is a common tool in statistical algorithms, but here we state it rather carefully 
since we later give a general method to convert any local search algorithm (that meets certain criteria) into 
one that finds all local maxima (see Section 4.2). 

Given that we can only ever hope for an additive approximation to a local maximum, one should be 
concerned about how the error accumulates when our goal is to find all local maxima. In fact, a naive 
strategy is to project onto the subspace orthogonal to the directions found so far, and continue in this 
subspace. However, such an approach seems to accumulate errors badly (the additive error of the last local 
maxima found is exponentially larger than the error of the first). Rather, the crux of our analysis is a novel 
method for bounding how much the error can accumulate (by refining old estimates) . 

Our strategy is to first find a local maximum in the orthogonal subspace, then run the local optimization 
algorithm again (in the original n-dimensional space) to "refine" the local maximum we have found. The 
intuition is that since we are already close to a particular local maxima, the local search algorithm cannot 
jump to some other local maxima (since this would entail going through a valley). 

4.1 Finding one local maximum 

Throughout this section, we will assume that we are given oracle access to a function f(u) and its gradient 
and Hessian. The procedure is also given a starting point u s , a search range fi, and a step size S. For 
simplicity in notation we define the following projection operator. 

Definition 4.1. Pio] ±u (v) = v — (u T v)u, Projj_ M (M) = M — (u T Mu)uu T . 

The basic step the algorithm is a modification of Newton's method to find a local improvement that makes 
progress so long as the current point u is far from a local maxima. Notice that if we add a small vector to u, 
we do not necessarily preserve the norm of u. In order to have control over how the norm of u changes, during 
local optimization step the algorithm projects the gradient V/ and Hessian IH(/) to the space perpendicular 
to u. There is also an additional correction term —d/d u f(u) ■ ||£|| 2 /2. This correction term is necessary 
because the new vector we obtain is (u + £)/ \\(u + £)|| 2 which is close to u — ■ u + ( + 0((3 3 ). Step 2 of 

the algorithm is just maximizing a quadratic function and can be solved exactly using Lagrangian Multiplier 
method. To increase efficiency it is also acceptable to perform an approximate maximization step by taking 
£ to be either aligned with the gradient Proj ±M V/(u) or the largest eigenvector of Projj_„(J{(/(u))). 

The algorithm is guaranteed to succeed in polynomial time when the function is Locally Improvable and 
Locally Approximable: 
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Algorithm 2. AllOPT, Input :/(«), /?, S, f3' , 5' Output: v u v 2 , v n , Vi \\vt - v*\\ < 7. 



1. Let vi = LoCALOPT(/,ei,/3,5) 

2. FOR i = 2 TO n DO 

3. Let gi be the projection of / to the orthogonal subspace of Vi, 1)2, ■ 

4. Let it' = LoCALOPT( 5 ,ei,/3',5'). 

5. Let vi = LocalOPT(/,u',/3,<5). 

6. END FOR 

7. Return ui, «2, v n 



Definition 4.2 ((7, /?, <5)-Locally Improvable). A function f(u) : R™ — > R is (7, /5, <5)-Locally Improvable, 
if for any u that is at least 7 far from any local maxima, there is a u' such that — u\\ 2 < /? and 

/(«')>/(«) + *■ 

Definition 4.3 <5)-Locally Approximablc). A function f(u) is locally approximablc, if its third order 
derivatives exist and for any u and any direction v, the third order derivative of / at point u in the direction 
of v is bounded by 0.01<5//3 3 . 

The analysis of the running time of the procedure comes from local Taylor expansion. When a function 
is Locally Approximable it is well approximated by the gradient and Hessian within a (3 neighborhood. The 
following theorem from [5] showed that the two properties above are enough to guarantee the success of a 
local search algorithm even when the function is only approximated. 

Theorem 4.4 ([5]). If \ f(u) — f*(u)\ < 5/8, the function f*(u) is (7, (3, 5) -Locally Improvable, f(u) is (/3,S) 
Locally Approximable, then Algorithm 1 will find a vector v that is 7 close to some local maximum. The 
running time is at most 0((n 2 +T) max/*/<5) where T is the time to evaluate the function f and its gradient 
and Hessian. 

4.2 Finding all local maxima 

Now we consider how to find all local maxima of a given function f*{u). The crucial condition that we need 
is that all local maxima are orthogonal (which is indeed true in our problem, and is morally true when using 
local search more generally in ICA). Note that this condition implies that there are at most n local maxima. 1 
In fact we will assume that there are exactly n local maxima. If we are given an exact oracle for /* and can 
compute exact local maxima then we can find all local maxima easily: find one local maximum, project the 
function into the orthogonal subspace, and continue to find more local maxima. 

Definition 4.5. The projection of a function / to a linear subspace S is a function on that subspace with 
value equal to /. More explicitly, if {vi,V2, ■■■,Vd} is an orthonormal basis of S, the projection of / to S is 
a function g : R d — > R such that g(w) = /(Xa=i w i v i)- 

The following theorem gives sufficient conditions under which the above algorithm finds all local maxima, 
making precise the intuition given at the beginning of this section. 

Theorem 4.6. Suppose the function f*(u) : R" — > R satisfies the following properties: 

1. Orthogonal Local Maxima: The function has n local maxima v* , and they are orthogonal to each other. 

2. Locally Improvable: f* is (7, j3, S) Locally Improvable. 

1 Technically, there are 2n local maxima since for each direction u that is a local maxima, so too is — u but this is an 
unimportant detail for our purposes. 
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3. Improvable Projection: The projection of the function to any subspace spanned by a subset of local 
maxima is (7', f3' , 8') Locally Improvable. The step size 8' > 10<5. 

4- Lipschitz: If two points \\u — u'\\ 2 < 3^/727, then the function value \f*{u) — f*(u')\ < <5'/20. 

5. Attraction Radius: Let Rad > 3-^/727+7', for any local maximum v* , let T be mm/*(u) for \\u — v*\\ 2 < 
Rad, then there exist a set U containing \\u — v*\\ 2 < 3-s/nr/ + 7' and does not contain any other local 
maxima, such that for every u that is not in U but is fj close to U , f*{u) < T . 

If we are given function f such that \f(u) — f*(u)\ < 8/8 and f is both ((3,6) and (ft 1 , 5') Locally 
Approximable, then Algorithm 2 can find all local maxima of f* within distance 7. 

To prove this theorem, we first notice the projection of the function / in Step 3 of the algorithm should 
be close to the projection of /* to the remaining local maxima. This is implied by Lipschitz condition 
and is formally shown in the following two lemmas. First we prove a "coupling" between the orthogonal 
complement of two close subspaces: 

Lemma 4.7. Given V\, V2, v k , each "/-close respectively to local maxima v*, v 2 , •••> v k (this is without 
loss of generality because we can permute the index of local maxima), then there is an orthonormal basis 
Wfe+i, Vk+2, v n for the orthogonal space of span{v\, V2, Vk} such that for any unit vector w G R" _fc , 
i= i WkVk+i *s 3V^7 close to 2^ i=1 w k v k+i . 

Proof: Let Si be span{vi,V2, S2 be span{v\ , v 2 , and S 1 ^, S 2 be their orthogonal subspaces 

respectively We first prove that for any unit vector v G Si , there is another unit vector v' G S 2 so that 
v T v' > 1 — Anj 2 . In fact, we can take v' to be the unit vector along the projection of v in S 2 . To bound the 
length of the projection, we instead bound the length of projection to S2. Since we know vfv' = for i < k 
and \\vi — < 7, it must be that (v*) T v' < 2j when 7 < 0.01. So the projection of v' in ^2 has length at 
most 2yjrvy and hence the projection of v' in S 2 has length at least 1 — 4wy 2 . 

Next, we prove that there is a pair of orthornormal basis {wfc+i, Ufc+2, v n } and {v*k+i, u*fc+2j •••) v* n } 
for Si and S 2 such that y?_~, w^Vk+i is close to X)"=i w k v *k+i- Once we have such a pair, we can 



simultaneously rotate the two basis so that the latter becomes v 

To get this set of basis we consider the projection operator to S 2 for vectors in S^- . The squared length 
of the projection is a quadratic form over the vectors in Si - So there is a symmetric PSD matrix M such 

that 



p roj g -L (v) 



= v T Mv for v G Si - Let {v k+ i,v k+2 , 



j} be the eigenvectors of this matrix M. As we 



showed the eigenvalues must be at least 1 — 8nj 2 . The basis for S 2 will just be unit vectors along directions 
of projections of Wj to S 2 . They must also be orthogonal because the projection operator is linear and 



n—k 



i — k 



WiPT0j s ±(v k+i 



i=l 



i — k 



J2\ 

i=l 



The second equality cannot hold if these vectors are not orthogonal. And for any w, 

T 

^ w k v* k +i = ^2 w 2 k (v k+i ) T v* k+i > 1 - 8n7 2 




So we conclude that the distance between these two vectors is at most 3y/nj. I 

Using this lemma we see that the projected function is close to the projection of /* in the span of the 
rest of local maxima: 

Lemma 4.8. Let g* be the projection of f* into the space spanned by the rest of local maxima, then \g*(w) — 
g(w)\ < 6/8 + 6' /20 < 6'/8. 
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Proof: The proof is straight forward because \g*(w) — g(w)\ < \f*{u) — f(u)\ + \f*(u) — f*(u')\ for some 
II u — 2 < 3-^/717, we know the first one is at most (5/8 and the second one is at most 6'/ 20 by Lipschitz 
Condition. ■ 

Now we are ready to prove the main theorem. 
Proof: [Theorem 4.6] By Theorem 4.4 the first column is indeed 7 close to a local maximum. We then 
prove by induction that if vi, V2, ■ Vk are 7 close to different local maxima, then Vk+i must be close to a 
new local maximum. 

By Lemma 4.8 we know gk+i is (7', /?',£') Locally Improvable, and because it is a projection of / its 
derivatives are also bounded so it is (/?',£') Locally Approximable. By Theorem 4.4 u' must be 7' close to 
local maximum for the projected function. Then since the projected space is close to the space spanned by 
the rest of local maxima, u' is in fact 7' + Sy/nj close to v% +1 (here again we are reindexing the local maxima 
wlog.). 

Now we use the Attraction Radius property, since u is currently in U , f*(u) > T, and each step we go 
to a point u' such that \\u' — u\\ < ft and f*(u') > f*(u) > T. The local search in Algorithm 1 can never go 
outside U, therefore it must find the local maximum u| , 1 . ■ 



5 Local search on the fourth order cumulant 

Next, we prove that the fourth order cumulant P*(u) satisfies the properties above. Then the algorithm 
given in the previous section will find all of the local maxima, which is the missing step in our main goal: 
learning a noisy linear transformation Ax + r\ with unknown Gaussian noise. We first use a theorem from 
[ : "'.] to show that properties for finding one local maxima is satisfied. 

Also, for notational convenience we set di = 2Da(uo)~ 2 and let d m in and d max denote the minimum 
and maximum values (bounds on these and their ratio follow from Claim 2.9). Using this notation P*{u) = 

Er=i*(« r fl?) 4 - 

Theorem 5.1 ([5]). When j3 < d mm / I0d max n 2 , the function P*(u) is (3^/3, (3, P* (u) /3 2 / 100) Locally 
Improvable and (/3, d m i n f3 2 /lOOn) Locally Approximable. Moreover, the local maxima of the function is 
exactly {±R*}. 

Proof: The proof appears in [5]. Here for completeness we show the proof using our notations. 

First we establish that P*{u) is Locally Improvable. Observe that this desirada is invariant under rotation, 
so we need only prove the theorem for P*(v) = Y^7=i^i v t- ^he gradient of the function is VP*(v) = 
A(divl,d 2 V2, ...,d n v^). The inner product of VP*(w) and v is exactly 4^" =1 divf = 4P*(v). Therefore the 
projected gradient <\> = Projj_„VP*(i;) has coordinate 4>i — 4:Vi(divf — P*(v)). Furthermore, the Hessian 
H = JC(P*(u)) is a diagonal matrix whose (i,i) th entry is \2diV 2 . 

Consider the case in which ||0|| > P*(v)(5/A. We can obtain an improvement to P*(w)/3 2 /100 because we 
can take £ in the direction of <fi and with ||£|| 9 = /3/20. The contribution of the Hessian term is nonnegative 
and the third term — 2P*(u) is small in comparison. 

Hence, we can assume ||^|| < P*(v)(3/4. Now let us write out the expression of ||0|| 2 

n 

Y,vf(d t v? - P*(v)) 2 < (3 2 (P*(v)) 2 /16. 

i=l 

In particular every term v 2 (divf — P*(v)) 2 must be at most f3 2 (P*(v)) 2 /16.. Thus for any i, either v 2 < (3 2 
or (d iV f - P*(v)) 2 < (P*{v)) 2 /16. 

If there are at least 2 coordinates k and / such that (diV 2 — P* (v)) 2 < (P* (v)) 2 /16, then we know for these 
two coordinates v 2 G [0.75P*(u)/dj, 1. 25 P*(v)/di}. We choose the vector £ so that = rv\ and £; = —rv k . 
Wlog assume £ • <j) > otherwise we use — £. Take r so that t 2 (v 2 + vf,) = [3 2 . Clearly ||£|| = (3 and £ • v = 

so £ is a valid solution. Also r 2 is lower bounded by fi 2 /(vf + v 2 .) > f p*{ u )(i%i+i/d k ) ' 
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Consider the function we are optimizing: 



<i> ■ i + l/2fM£ - 2P*(u) ||£|| a > 1/2£ T £T£ - 2P*(u)(3 2 = 6r 2 v 2 k vf(d k + d,) - 2P»/3 2 

- f r2p * (u)2 ^T i " 2P * (M)/?2 - ^ p »^ 

In the remaining case, all of the coordinates except for at most one satisfy vf < 1 . Since we assumed 
/3 2 < -, there must be one of the coordinate v k that is large, and it is at least 1 — nf3 2 . Thus the distance 
of this vector to the local maxima is at most 3y/n(3. I 

We then observe that given enough samples, the empirical mean P'{u) is close to P* (it). For concentration 
we require every degree four term ZiZjZ k zi has variance at most Z . 

Claim 5.2. Z = 0(d 2 min X min (A) s ||S||| + d 2 nm ). 

Proof: Wc will start by bounding E[(z t ZjZ k zi) 2 } < ~E[(zf+z®+zl+zf)}. Furthermore E [zf] < 0(E[(B~ 1 Ax) 1 * + 
(_B _1 r;)f]). Next we bound ~E[(B~ 1 i])f], which is just the eighth moment of a Gaussian with variance at most 
|| J B- 1 S J B- T || 2 < WB^WlW^h < d^ n A min (A)- 2 ||S|| 2 . Hence we can bound this term byO(\\B- l EB- T \\i) = 
0(d 2 , lin A m i„(^4) 8 |jS|||). Finally the remaining term E[(i3 _1 Ar)f ] can be bounded by 0{d 2 nin ) because the 
variance of this random variable is only larger if we instead replace x by an n-dimensional standard Gaussian. 
■ 

Lemma 5.3. Given 2N samples j/i, y2, Uni Vi-, •••) d'n> su PP ose columns of R' = B^ 1 ADa{uq) 1 ^ 2 are 
€ close to the corresponding columns of R* , with high probability the function P'(u) is 0(d max n 1 ^ 2 e + 
n 2 (iV/Zlogn) -1 ' 2 ) close to the true function P*(u). 

Proof: P'{u) is the empirical mean of F(u,y,y') = —(u T B~ 1 y) A + 3(it T B~ 1 y) 2 (u T B~ x y') 2 . In Section 2 
wc proved that P'(u) = E y/y < F(u, y, y') = Yl?=i (u T Ri) 4 = S"=i ^i(u T Ri) 4 - First, we demonstrate 

that P'(u) is close to P*(u), and then using concentration bounds we show that P'(u) is close to P'(u) (with 
high probability) over all u. 

The first part is a simple application of Cauchy-Schwartz: 



\P'(u) - P*(u)\ = [(u T R>) - (u T R*)] ■ [(u r i$ + u T Rl){{u T R! i ) 2 + (u T R;) 2 )] 

i=l 



Y,^ T {R' t - R*)) 2 ■ (3 \\u T R' + u T R* || 2 ) < 6d max n 1/2 e. 

i=l 

The first inequality uses the fact that ((u T i?£) 2 + (u T R*) 2 ) < 3, the second inequality uses the fact that 
when e is small enough, ||u T i?'|| 2 < 2. 

Next we prove that the empirical mean P'(u) is close to P'(u). The key point here is we need to prove 
this for all points u since a priori wc have no control over which directions local search will choose to explore. 
Wc accomplish this by considering P'(u) as a degree-4 polynomial over u and prove that the coefficient of 
each monomial in P'(u) is close to the corresponding coefficient in P'(u). This is easy: the expectation of 
each coefficient of F(u,y,y') is equal to the correct coefficient, and the variance is bounded by O(Z). The 
coefficients are also sub-Gaussian so by Bernstein's inequality the probability that any coefficient of P'(u) 
deviates by more than e' (from its expectation) is at most e~ n ^ e N / z \ Hence when TV > 0{Z logn/e' 2 ) with 
high probability all the coefficients of P'(u) and P'(u) are e' close. So for any u: 

n 

\P'(u)-P'(u)\<e'{Y.W l \) A <z'n 2 . 

i=l 



\ 
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Algorithm 3. Recover, Input:P, P'(u), R, e Output: A, E 

1. Let Da(u) be a diagonal matrix whose i entry is | (P'(-Ri)j 

2. Let A = BRDAiu) 1 ^ 2 - 

3. Estimate C = E[yy T ] by taking 0((||A|| 2 + ||E|| 2 ) 4 n 2 e" 2 ) samples and let 6 = ^ E^Li 2/«yf- 

4. Let £ = C - A4 T 

5. Return .A, £ 



Therefore P'(u) and P*(u) are 0(d max n^ 2 e + n 2 (N/Z log n)~^ 2 ) close. ■ 

This proof can also be used to show that the derivatives of the function P'(u) is concentrated to the 
derivatives of the true function P*(u) because the derivatives are only related to coefficients, therefore P'(u) 
is also (j3, d m i n 1 jlQQn) Locally Approximable. 

The other properties required by Theorem 4.6 are also satisfied: 

Lemma 5.4. For any \\u — u'\\2 < r, \P*(u) — P*(u')\ < ^dmax^^r. All local maxima of P* has attraction 
radius Rad > d min / I00d max . 

Proof: The Lipschitz condition follows from the same Cauchy-Schwartz as appeared above. When two 
points u and u' are of distance r, \P*(u) — P*(u')\ < hd max n x l 2 r . Finally for the Attraction Radius, we 
know when Zyfnrf + 7' < d m i n /100d max , we can just take the set U to be u T R* > 1 — d min /50d max . For 
all u such that u T R* G [1 — d min /25d max , 1 — d min /50d, nax ] (which contains the /3 neighborhood of U), we 
know the value of P*(u) < T. U 

Applying Theorem 4.6 we obtain the following Lemma (the parameters are chosen so that all properties 
required are satisfied): 

Lemma 5.5. Let ft = Q{{d min /d max ) 2 ), (3 = mm{'ynr 1 ^ 2 ,n((d min /d m ax)' i n~ 3 - 5 )}, then the procedure 
Recover(/, /3, d m i n /3 2 /100n , /?', d m i n (3' 2 /100n) finds vectors V\,V2, ■■■■,v n , so that there is a permutation 
matrix II and ki G {±1} and for all i: \\vi — (RHDiag(ki))*\\ 2 < 7. 

After obtaining R = [v\, V2, v n ] we can use Algorithm 3 to find A and E: 

Theorem 5.6. Given a matrix R such that there is permutation matrix II and ki G {±1} with \\Ri — 
ki(R*TL)i\\2 < 7 for all i, Algorithm 3 returns matrix A such that \\A— ATLDiag{ki)\\F < 0(7 || A\\ 2 n 3 ! 2 j A m j n (A)). 
Ifl < 0{e/ \\A\\ln 3 ' 2 \ mm (A)) x min{l/ ||A|| 2 , 1}, we also have ||E-E|| F < e. 

Recall that the diagonal matrix Da(u) is unknown (since it depends on A), but if we are given R* (or 
an approximation) and since P*(u) = ^2™ = idi(u T R*) 4 , we can recover the matrix Da(u) approximately 
from computing P*(i?*). Then given Da(u), we can recover A and E and this completes the analysis of our 
algorithm. 

Proof: By Lemma 2.11 we know the columns of R' is close the the columns of R (the parameters will 
be set so that the error is much smaller than 7), thus \\Ri — ki(R'Ii)i\\2 < 7. Applying Lemma 5.3 we 
obtain: \P'(R 1 ) — P*(P;)| <C 7- Furthermore, when \\Ri — hR^^,^^ < 7 we know that P*(Pi)/d n -i(,:) G 
[1 — 37, 1 + 37] (here we are abusing notation and use the permutation matrix as a permutation). Hence 
A4 («)»,*/ (D A (u)) n -i {l) n -i {z) G [1 - 37, 1 + 37]. We have: 

A\ = BR t D A {u)ll /2 and (ALTDiagCfc;)), = BR^ 1(i) (JM«))n-i(<),n-*(0 
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and their difference is at most 0(-/ \\B\\ 2 (£) > i(u)) n _'i^ n -i^\)- Hence we can bound the total error by 
0(7 115112 HDa^)- 1 / 2 !^). We also know \\B\\ 2 < \\A\\ 2 WD^u) 1 / 2 ^ because BB T AD A (u)A T , so this 
can be bounded by 0(7 |[A|| 2 HD^ii)^ \\Da{u)~ 1 ^ 2 \\f)- Applying Claim 2.9, we conclude that (with high 
probability) the ratio of the largest to smallest diagonal entry of Da(u) is at most n 2 \\A\\ 2 2 / X m i n (A) 2 . So 
we can bound the error by 0(7 ||^4||2 n 3 ^ 2 /X m i n (A)). 

Consider the error for S: Using concentration bounds similar but much simpler than those used in 
Lemma 5.3, we obtain that j|C - C\\ F < l/2e, so ||S - H\\ F < \\C - C\\ F - \\AA T - AA T \\ F < e/2 + 
2 \\A\\ 2 ||AnDiag(fc,) -A\\ F + ||AnDiag(fc l ) - A\\% < e. ■ 

Conclusions 

ICA is a vast field with many successful techniques. Most rely on heuristic nonlinear optimization. An 
exciting question is: can we give a rigorous analysis of those techniques as well, just as we did for local 
search on cumulants? A rigorous analysis of deep learning — say, an algorithm that provably learns the 
parameters of an RBM — is another problem that is wide open, and a plausible special case involves subtle 
variations on the problem we considered here. 
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